Chapter 7 Diversity

load("data/data.Rdata")

7.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  filter(genome %in% labels(dist)[[1]]) %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == sample)) %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Species",
          breaks=c("Pk","Eb","Ha"),
          labels=c("Pipistrellus kuhlii","Eptesicus bottaew","Hipsugo ariel"),
          values=c("#e5bd5b", "#6b7398","#e2815a")) +
      scale_fill_manual(name="Species",
          breaks=c("Pk","Eb","Ha"),
          labels=c("Pipistrellus kuhlii","Eptesicus bottaew","Hipsugo ariel"),
          values=c("#e5bd5b50", "#6b739850","#e2815a50")) +
      facet_wrap(. ~ metric, scales = "free", ncol=4) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank())

7.2 Beta diversity

beta_q0n <- genome_counts %>%
  select(where(~!all(. == 0))) %>% # remove empty samples
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  select(where(~!all(. == 0))) %>% # remove empty samples
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts %>%
  select(where(~!all(. == 0))) %>% # remove empty samples
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts %>%
  select(where(~!all(. == 0))) %>% # remove empty samples
  filter(genome %in% labels(dist)[[1]]) %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

7.2.1 Richness (q0n)

tinytable_5x9jjr762wsc044wvwhz
Homogeneity of variances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
2 0.0005478071 0.0002739036 0.3643863 999 0.695
69 0.0518662354 0.0007516846 NA NA NA
tinytable_bsqz9umml87do4s6lx1h
Permanova
term df SumOfSqs R2 statistic p.value
species 2 0.004372878 0.03552136 1.270621 0.239
Residual 69 0.118732719 0.96447864 NA NA
Total 71 0.123105597 1.00000000 NA NA
beta_q0n$C %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
    scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

7.2.2 Neutral (q1n)

tinytable_bvme5sif0eno4uqxxsa5
Homogeneity of variances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
2 0.2034005 0.10170027 7.084138 999 0.003
69 0.9905677 0.01435605 NA NA NA
tinytable_0dea5y7c62g6d80b4t8a
Permanova
term df SumOfSqs R2 statistic p.value
species 2 1.567475 0.06477038 2.389336 0.003
Residual 69 22.633010 0.93522962 NA NA
Total 71 24.200485 1.00000000 NA NA
beta_q1n$C %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
    scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

7.2.3 Phylogenetic (q1p)

tinytable_7dqt55d6chl9khxvsvnl
Homogeneity of variances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
2 0.135256 0.06762802 3.065569 999 0.057
69 1.522175 0.02206051 NA NA NA
tinytable_8w678d19vg046w0zbgus
Permanova
term df SumOfSqs R2 statistic p.value
species 2 0.4397545 0.04243533 1.528898 0.139
Residual 69 9.9231780 0.95756467 NA NA
Total 71 10.3629324 1.00000000 NA NA
beta_q1p$C %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
    scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

7.2.4 Functional (q1f)

tinytable_ffv2i50npvtyircyt56o
Homogeneity of variances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
2 0.2176262 0.10881312 2.510851 999 0.098
69 2.9902627 0.04333714 NA NA NA
tinytable_qd2geuz0c8si0pd1rvkp
Permanova
term df SumOfSqs R2 statistic p.value
species 2 0.9379128 0.0763021 2.849874 0.044
Residual 69 11.3541841 0.9236979 NA NA
Total 71 12.2920969 1.0000000 NA NA
beta_q1f$C %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace = 0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(species) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = species)) +
    scale_color_manual(values = c("#e5bd5b", "#6b7398","#e2815a")) +
    scale_shape_manual(values = 1:10) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )